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Study of wave chaos in a randomly-inhomogeneous oceanic acoustic waveguide: 
spectral analysis of the finite-range evolution operator 
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The proplem of sound propagation in an oceanic waveguide is considered. Scattering on random 
inhomogeneity of the waveguide leads to wave chaos. Chaos reveals itself in spectral properties of 
the finite-range evolution operator (FREO). FREO describes transformation of a wavefield in course 
of propagation along a finite segment of a waveguide. We study transition to chaos by tracking 
variations in spectral statistics with increasing length of the segment. Analysis of the FREO is 
accompanied with ray calculations using the one-step Poincare map which is the classical counterpart 
of the FREO. Underwater sound channel in the Sea of Japan is taken for an example. Several 
methods of spectral analysis are utilized. In particular, we approximate level spacing statistics 
by means of the Berry-Robnik and Brody distributions, explore the spectrum using the procedure 
elaborated by A. Relano with coworkers (Relano et al, Phys. Rev. Lett., 2002; Relano, Phys. 
Rev. Lett., 2008), and analyze modal expansions of the eigenfunctions. We show that the analysis 
of FREO eigenfunctions is more informative than the analysis of eigenvalue statistics. It is found 
that near-axial sound propagation in the Sea of Japan preserves stability even over distances of 
hundreds kilometers. This phenomenon is associated with the presence of a shearless torus in the 
classical phase space. Increasing of acoustic wavelength degrades scattering, resulting in recovery 
of localization near periodic orbits of the one-step Poincare map. Relying upon the formal analogy 
between wave and quantum chaos, we suggest that the concept of FREO, supported by classical 
calculations via the one-step Poincare map, can be efficiently applied for studying chaos-induced 
decoherence in quantum systems. 

PACS numbers: 05.45.Mt, 43.30.Cq, 03.65.Yz, 05. 45. Ac, 43.30.Ft 



I. INTRODUCTION 

Sound speed in the deep ocean typically has a mini- 
mum at some depth. This results in formation of a re- 
fractive waveguide, the so-called underwater sound chan- 
nel, which prevents sound waves from contact with the 
absorbing bottom. As sound absorption within water col- 
umn is fairly weak, an underwater sound channel enables 
sound propagation over distances of thousands kilome- 
ters. The largest distance had been achieved using ex- 
plosive charges in the seminal experiment on sound trans- 
mission from Perth to Bermuda in 1960 [l|, [Ji- 
lt is realized that small sound-speed variations induced 
by oceanic internal waves lead to Lyapunov instability 
and chaos of sound rays. In the mathematical sense, ray 
chaos is an analogue of classical chaos in Hamiltonian 
systems. Indeed, ray motion in a waveguide is equiv- 
alent to motion of a point particle in a potential well, 
and sound-speed variations along a waveguide play the 
role of a nonstationary perturbation. Reciprocal Lya- 
punov exponent for chaotic rays typically is of about 
several tens kilometers \3\ , therefore, the problem of ray 
chaos is mainly important for long-range sound propa- 
gation. During the last two decades ray chaos in ocean 
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acoustics was an object of intense research, both theo- 
retical and experimental I^tSI- Considerable attention 
was paid to wave chaos [3|, I9l4l4| . The term wave chaos 
relates to wavefield manifestations of ray chaos. It was 
found that interference makes wave refraction more regu- 
lar than it is anticipated from ray modeling, albeit influ- 
ence of ray chaos persist even for very low sound frequen- 
cies |15| . This problem becomes especially important due 
to the growing interest to hydroacoustical tomography, 
i.e. monitoring of environment using sound signals. The 
classical scheme of tomography developed by Munk and 
Wunsch [l6| is based on computation of eigenrays con- 
necting the source and the receiver. It was shown in 
[17t that ray chaos leads to exponential proliferation of 
eigenrays with increasing distance. As a result, the in- 
verse problem becomes ill-posed, impeding environment 
reconstruction. However, wave-based corrections "sta- 
bilize" wave refraction, i. e. the standard semiclassical 
approximation typically overestimates ray chaos. Thus, 
one needs either an improved version of the semiclassical 
approximation for proper computation of eigenrays, or 
some approach for making implications about eigenray 
stability relying upon wave modeling, whereby a priori 
taking into account the wave-based suppression of ray 
chaos. 

Theory of ray and wave chaos extensively exploits 
methods borrowed from the theory of Hamiltonian dy- 
namical systems [SJ, [l8|. In particular, ray motion is 



often studied by means of the phase space representa- 
tion. This provides the clear geometric interpretation 
of ray dynamics for toy models of deterministic range- 
periodic waveguides. For instance, one can easily sepa- 
rate domains of initial conditions corresponding to sta- 
ble and chaotic rays from each other. However, realistic 
underwater sound channels are not range-periodic, and 
their range inhomogeneity should be rather described 
as a stochastic process. On the other hand, statistical 
methods [T^ [20| are implicitly based on the assump- 
tion of ergodic chaos and ignore the existence of phase 
space domains of finite-range stability [2l|, ^^. There- 
fore, it is reasonable to elaborate some theoretical ap- 
proach combining advantages of statistical and determin- 
istic methods. Recently it had been shown that a bridge 
between the deterministic and statistical descriptions of 
wave propagation in random media can be built up by 
introducing the finite-range evolution operator (hereafter 
FREO) [14, 23]. FREO acts as a propagator determining 
wavefield transformation between two vertical sections of 
a waveguide. Spectral analysis of the FREO in terms of 
the random matrix theory allows one to estimate contri- 
bution of ray chaos to wave dynamics and track transi- 
tion to chaos with increasing distance between the sec- 
tions. A somewhat different approach had been offered 
in [24l|. These approaches allow one to generalize the 
well-developed spectral theory of quantum chaos (see, for 
instance, |25| ) on one-way wave propagation in random 
media. On the other hand, concept of the FREO can 
serve as a promising method for studying chaos-induced 
decoherence in randomly-driven quantum systems |26| . 

In the present paper we study spectral statistics of the 
FREO for the Sea of Japan. Attention is concentrated 
on the track connecting the Gamov peninsula and Kita- 
Yamato bank. The length of the track is about 350 km. 
Our interest to this waveguide is motivated by results 
of the experiment conducted there in 2006 \2T\ , indicat- 
ing on high stability of near-axial propagation. Simi- 
lar behavior was observed in a earlier experiment with 
a slightly different propagation track z8]. It should be 
noted that stability of near-axial propagation is atypical 
for the deep ocean [J, [Ij] . 

The paper is organized as follows. The next section 
represents basic equations describing long-range sound 
propagation in the ocean. In Section Hill we describe the 
waveguide used in the paper. Section lTVl contains descrip- 
tion of the FREO and methods of its analysis. Section 
rvlis devoted to the classical counterpart of the FREO, 
namely the one-step Poincare map. In Section I VII we 
perform statistical analysis of the FREO for the under- 
water sound channel in the Sea of Japan. In Conclusion 
we summarize and discuss the results obtained. 



II. GENERAL EQUATIONS 



to reduce the initial three-dimensional problem of wave 
propagation in an underwater sound channel to a two- 
dimensional one by assuming cylindrical symmetry and 
neglecting azimuthal coupling. Sound refraction is gov- 
erned by spatial variability of sound speed 



c{z, r) = co + Ac{z) + Sc{z, r), 



(1) 



where z is depth, r is range coordinate, cq is a reference 
sound speed. Sound-speed variations obey the double 
inequality 



l-^cl: 



< lAcI, 



<co. 



(2) 



Left inequality implies that the range-dependent term 
can be treated as a weak perturbation of the background 
sound-speed profile. This term is mainly contributed 
from oceanic internal waves. Right inequality means that 
variations of the refractive index are weak, and only those 
waves which propagate with small angles with respect to 
the horizontal plane can avoid contact with the absorbing 
ocean bottom. Thus, one can invoke the small-angle ap- 
proximation, in which an acoustic wavefield is governed 
by the standard parabolic equation 



i 9$ 
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+ [U{z) + eV{z,r)]<^, (3) 



where wave function $ is linked to acoustic pressure u 
by means of the formula u — $exp(ifco'')/v^- Here the 
denominator ^/r responds for the cylindrical spreading of 
sound. Quantity fco is the reference wavenumber related 
to sound frequency f as kQ — 27r//co. Functions U{z) 
and V{z, r) are determined by spatial sound-speed varia- 
tions. In the the small-angle approximation they can be 
expressed as 

rrr N Ac(z) Scjz, r) 

U{z) = , V(z, r) = . (4) 
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According to ([2]), |y|max < |C/|max, that is, function 
V{z, r) can be treated as a small perturbation. One can 
easily see that the replacement 



kn^ ^ h. 



t 



(5) 



transforms the parabolic equation ([3|) into the 
Schrodinger equation for a particle with unit mass. 
This circumstance enables study of wave propagation 
using the approaches developed in quantum mechanics. 
In this relationship, function U{z) serves as an unper- 
turbed potential. As r is a timelike variable, V{z, r) 
plays the role of a nonstationary perturbation. 

In the short-wavelength limit kg ^ 00, solution of the 
parabolic equation ([3]) can be expressed as a sum of rays 
whose trajectories are governed by the Hamiltonian 



Ocean is a layered media, and its horizontal variabil- 
ity is much weaker than vertical one. This allows one 
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p^ 



+ U{z) + Viz,r), 



(6) 



p, 



dp _ dH _ dU dV 
dr dz dz dz 



where p = tan^, X is ray grazing angle (i. e. angle with 
respect to the horizontal plane). The respective Hamil- 
tonian equations read 

dz _ dH 
dr dp 

Owing to the analogy with classical mechanics, p is re- 
ferred to as ray momentum. 

III. MODEL OF A WAVEGUIDE 

A. Background sound-speed profile 
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FIG. 1. (a) Unperturbed sound-speed profile, (b) the first 
empirical orthogonal function of the sound-speed perturba- 
tion (solid) and its smoothed approximation (dotted). 



Model of the underwater sound channel in the Sea of 
Japan was elaborated using the hydrological data from 
the database [231 ■ Function U{z) corresponding to the 
background sound-speed profile was approximated by the 
expression 



U{z) = 



Ui{z), z ^ zo, 

[/2(z), Z > Zq. 



(8) 



[7^(^)^£le-M, 
cq 
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(9) 



Co = 1455 m/s, ci = 70 m/s, zq — 250 m is the depth of 
the channel axis, i. e. the depth with the minimal sound 
speed, zi = 30 m, g = 0.017 s^^. (see Fig. [J). The ocean 
bottom is assumed to be flat and placed at the depth 
h = 3 km. We consider only the deep-water propaga- 
tion, albeit the source in the aforementioned experiments 
[271 128| was mounted into the bottom in the coastal zone 
near the Gamov peninsula. The shallow-water part of 
the waveguide was relatively short, less than 30 km, and 
didn't have significant bottom features which could re- 
markably alter ray stability. 
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FIG. 2. Ray cycle length vs initial ray momentum for the 
source located at the channel axis. 

Expressions ([8]) and ([9|) permit analytical derivation of 
the basic model characteritics in the absence of horizontal 
inhomogeneity. One can introduce the ray action 



1 



/=- / ^2[E-Uiz)]dz, 



(10) 



where Zmin and z^ax are the upper and lower ray turning 
points, respectively, and E = 0.5p^ + U{z). Ray action 
describes steepness of a ray trajectory and enters into the 
Einstein-Brillouin-Keller quantization rule. 



koI„ 



1/2, 



1,2,. 



(11) 



establishing the link between normal modes of the unper- 
turbed waveguide and modal rays. Here Im is the action 
of a modal ray, and both ray turning points are assumed 
to be inside the water column, that is, rays undergo total 
internal reflection due to smooth vertical gradient of the 
refractive index n = co/c(z, r). Normal modes are the 
solutions of the Sturm-Liouville equation 



1 d^njz) 
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+ U{z)(t)n{z) = En(t>n{z), 



(12) 



with appropriate boundary conditions at the ocean sur- 
face and bottom. The main contribution to the mode m 
is given from the rays with / ~ /,„ [9]. Integration of 
([TU| yields 



I = 2ziy/2E~ V^hi(V?+V£-l)-Ve-l 
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where 



E/Er, 



En 



(ci/co)exp(-zo/2;i). Ray 



cycle length, i. e. horizontal distance between two suc- 
cessive upper (or lower) ray turning points, can be deter- 
mined as 



D(E) = 27r— 
^ ' dE 



2zl^/|ln(Vi+^/F^T) + ?^H^/^^VI^T. (14) 
V -c/ g 

Fig. [2] represents dependence of ray cycle length on ini- 
tial ray momentum for the source located at the channel 
axis z = 250 m. Initial momentum po depends on E 
as po — \/2(E — i?min)- Function D{pq) is nonmono- 
tonic and has two extrema, the sharp maximum and the 
smooth minimum. The latter one can give rise to a so- 
called weakly-divergent beam [30l - t33| . Its low divergence 
is associated with approximate equality of cycle length 
values for rays forming the beam. It will be demon- 
strated in Section |V] that the local minimum of D{po) 
plays a significant role in ray stabilty. 



B. Sound-speed perturbation 

Model of the internal-wave-induced sound-speed per- 
turbation was built up in several steps. Firstly, we com- 
puted the range-averaged profile of buoyancy frequency, 
using the hydrological data [2g|. Then, we calculated 
realizations of the sound-speed perturbation using the 
method proposed in J34i] . In order to facilitate numerical 
simulation, the perturbation was expanded over empiri- 
cal orthogonal functions [3a| 

Sciz, r) = {dc{z)) + ^ b„{r)Y„{z). (15) 

n 

Empirical orthogonal functions Yn{z) are the eigenvec- 
tors of the covariance matrix K with elements 



K. 



1 ^ 

-J2[Sc,{z,) ~ {Sc{z,)Wci{zj) - {Sc{z,))], (16) 



1=1 



where index I numbers L statistically independent real- 
izations of dc{z), {zi} is a vector of depth values, and an- 
gular brackets mean ensemble average. As 5c is caused by 
internal waves, one can set (Sc) = 0. Eigenvalues of the 
matrix K quantify contributions from the corresponding 



eigenvectors in the expansion ((T5|) . It was found that the 
contribution of the first orthogonal function prevails, and 
one can fairly represent the sound-speed perturbation as 
the product 



V{z,r)^h{r)Yi{z) 



(17) 



where bi{r) is a random function. For simplicity, it 
is assumed that 61 (r) is a stochastic process with the 
exponentially-decaying autocorrelation function 



(6i(r)6i(r'))=exp(-|r-r'|/f), 



(18) 



where the correlation length f is taken of 10 km, that is 
typical for the deep ocean [3^. Then the realizations of 
61 (r) can be computed via the formula &i — \/2fr], where 
77 is a solution of the Ornstein-Uhlenbeck stochastic dif- 
ferential equation 



dr] 1 . N 1 ,/ N 

dr r r 



(19) 



The procedure of solving this equation is described, for 
instance, in [37| . Here ^ is a Gaussian white noise satis- 
fying 



(e(r))=0, (e(r)e(r'))='5(r-r'). 



(20) 



The resulting function fei(r) satisfies {h\{r)) ~ 1. 

The calculated function Yi{z) is depicted in Fig.[TJb). 
It involves step-like changes in the depth interval from 
100 to 300 meters. These changes are caused by the 
depth discretization of the hydrological data and phys- 
ically irrelevant. Wave modeling is insensitive to them, 
but ray-based calculations can be significantly affected. 
Therefore, in ray calculations we use a smooth analytical 
approximation 



Yi ^ Ay^eM~v2) + Sexp(-y^), 



(21) 



where A = 0.0027, j/a — z/za., n 
B = 2- 10~^, y^ = {z- Zc)/z\„ Zb = 



= 1.1, Za = 24 m, 
50 m, Zc = 200 m. 



IV. FINITE-RANGE EVOLUTION OPERATOR 

Finite-range evolution operator (FREO) had been in- 
troduced for studying wave propagation in a randomly- 
inhomogeneous waveguide in [IJ, |23|. Its quantum- 
mechanical analogue was earlier considered in }2u\. Ba- 
sically, a finite-range evolution operator (FREO) is an 
element of one-parameter group generated by the opera- 
tor in the right-hand side of the equation Q. Consider 
a solution $(z, r) of the parabolic equation ^ comple- 
mented with the standard boundary conditions 



*i 



az 



= 



(22) 



z^h 



and the initial condition $(2;, r 
belongs to i^[0, h] and satisfies 



0) 



= $(z), where $(z) 
Then FREO G{t) 



is defined on the subspace of L^[0, h] (restricted by (|22)) ) 
as 



G{t)^{z) = $(z, 



(23) 



By definition, the FREO describes transformation of a 
wavefield in course of propagation along a finite waveg- 
uide segment of length r. Each realization of inhomo- 
geneity produces its own realization of the FREO. Our 
interest is concerned with statistical properties of the 
FREO and their connection to classical ray stability. 

Note that the choice of the hard wall boundary con- 
dition at the bottom (|22p is typical for the deep ocean 
acoustics problems when the attention is restricted to the 
trapped modes, whose propagation is not affected by the 
bottom interaction. Under these conditions no energy is 
absorbed by the bottom. Hence, if the attenuation in the 
sea water is negligible (this is true for the sound frequen- 
cies of our interest) and refraction index in ([3]) has no 
imaginary part, then the FREO is a unitary operator. 

FREO can be represented as a matrix in the basis of 
normal modes 4>j{z) satisfying the Sturm-Liouville prob- 
lem (|12I22|) . Matrix elements of the FREO are given by 



Gmnir) = / 0™(0)G(r)0„(z)dz. (24) 



Thus, the matrix elements Gmn are complex-valued am- 
plitudes of modal transitions. For the range-independent 
waveguide, the matrix of FREO is diagonal with 

Eigenvalues and eigenvectors of the FREO obey the 
equation 



G*m(z, r) = gm'^miz, r). 



(25) 



Owing to unitarity, eigenvalues gm can be recast as 

ff™ = e-*'=°^'", e^GSft. (26) 

Since eigenvalues of the FREO belong to the unit circle in 
the compl ex p lane, the FREO corresponds to the circular 
ensemble [25| . The FREO has much in common with the 
Floquet operator governing wave propagation in a range- 
periodic waveguide [111 |l2| and quantum dynamics in 
time-periodic systems. For instance, quantity £,„ is the 
analogue of quasienergy in quantum mechanics. 

Note that eigenvalues of G{t) may be easily computed 
using its matrix representation Gmn{T). To accomplish 
this one has to clip a finite block of this (infinite) matrix 
corresponding to the trapped modes, neglecting their in- 
teraction with high-order modes. This simplification is 
reasonable and does not affect accuracy of the eigenval- 
ues computation (and the numerics confirms that) since 
the prevailing small-angle propagation corresponds to the 
low-order modes. 



A. Level spacing statistics 

Wave chaos reveals itself in the spectrum of the FREO, 
in particular, in the statistics of level spacings. A level 
spacing is defined as 



fcoM(e, 
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(27) 



where tm increases with increasing rra. A/ is the to- 
tal number of eigenvalues for a single realization of the 
FREO, equal to the number of trapped modes. 

Level spacing statistics can be studied in terms of the 
random matrix theory |25|. Regular dynamics implies 
that the matrix of the FREO consists of separate inde- 
pendent blocks. Then level sequences contributed from 
different blocks are statistically independent, therefore, 
the resulting level spacing distribution obeys the Poisson 
law 



p(s) ~ exp(-s) 



(28) 



Under conditions of ergodic chaos, all normal modes are 
coupled to each other. This results in repulsion of neigh- 
bouring levels, the phenomenon closely related to spec- 
tral splittings induced by tunneling [38|. In this case 
level spacing statistics is described by the Wigner-Dyson 
distribution 



P(s) 



s^ exp I 



-Gs' 



(29) 



where constants C and G depend on symmetries of the 
FREO. As the FREO doesn't possess the symmetry 
r — > — r, it corresponds to the circular unitary ensem- 
ble (CUE) with C = 2 and C = 4/7r [H]. 




FIG. 3. Berry-Robnik distribution with v-c — 0.9, Vi — 0.5 
and Vr ~ 0.1. 



The most interesting case is mixed phase space, with 
the coexistence of regular and chaotic domains. Then 



level spacing statistics should be described by some com- 
bination of Poisson and Wigner-Dyson laws. In the short- 
wavelength limit one can use the Berry-Robnik distribu- 
tion (si 



Pis) 



' erfc 
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v^ s I exp 
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2 2 



exp(— Wrs), (30) 



where v^ and Vc are relative phase space volumes corre- 
sponding to regular and chaotic ray motion, respectively, 
Vj: + Vc = I. In (j30p . it is assumed that phase space con- 
sists of only two distinct domains: one regular and one 
chaotic. 

Taking into account the correspondence between rays 
and modes established by the WKB representation ([TTj). 
one can regard v^- as a fraction of modes whose refraction 
is regular, and Vc as a fraction of modes exhibiting wave 
chaos. In the limiting cases v^ = 1 and Wc = 1, Berry- 
Robnik formula (pO|) reduces to the Poisson distribution 
and the Wigner-Dyson distribution for the orthogonal 
ensemble (C = 1), respectively. As is shown in Fig. [31 
Berry-Robnik distribution undergoes smooth transition 
from the Poisson to the Wigner-Dyson law as v^ decreases 
from 1 to 0. Wigner-Dyson distributions for orthogonal 
{( = 1) and unitary (C = 2) ensembles are close to each 
other, therefore, one can use (j30p as an approximation 
for the circular unitary ensemble. Berry-Robnik distri- 
bution is based on the assumption of the total statisti- 
cal independence of the matrix blocks corresponding to 
regular and chaotic dynamics. This assumption is com- 
pletely fulfilled only in the semiclassical limit. As wave 
corrections grow, independence degrades due to regular- 
to-chaotic tunneling [40|. Hence, Berry-Robnik formula 
cannot work perfectly for low-frequency sound propaga- 
tion. In this case, a better fit is expected with the Brody 
distribution [41i] 



p(s) = (/3-t-l)A/3S^exp(-A^s^+i), 



(31) 



where A^ — [r(4xx)]'^^^; ^ ^^ ^^^ Euler gamma function. 
(3 — corresponds to the Poisson distribution, /3 = 1 
yields the Wigner-Dyson distribution with C = 1. Un- 
fortunately, the Brody distribution is semiempirical and 
doesn't have an explicit physical interpretation in the in- 
termediate regime < /3 < 1. 

Functions ([50)1 and ([51]) should describe level spacing 
statistics for single realizations of the FREO. Very un- 
fortunately, this is problematic because one encounters 
insufficiency of the statistical ensemble. Indeed, long- 
range sound propagation is feasible only with low acous- 
tic frequencies of tens or hundreds Hz. In this frequency 
range, number of trapped modes doesn't exceed several 
hundreds. To resolve this problem, we consider ensemble- 
averaged level spacing distribution 



pis, t) 
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n=l 
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(32) 



where Pnis^ t) is a level spacing distribution correspond- 
ing to n-th realization of FREO. Fitting the function 
p{s, t) with the Berry-Robnik distribution (pO| . one can 
estimate number regularly propagating modes for var- 
ious values of t and, whereby, track the transtion to 
chaos with increasing r. However, it should be noted 
that formula ([32]) enables accurate estimate of Vi- only 
if fluctuations of Wj- are weak, otherwise one should take 
into account nonlinearity of p as a function of Vi- in (|30p . 
Likewise we can observe transformation of level spacing 
statistics by fitting (15^ with the Brody distribution (|31l) . 
Then the value of (3 corresponding to the best fit should 
increase with increasing r from to 1, reflecting wave- 
field stochastization. Such calculations were performed 
in [14l.l23l| for a waveguide with a perturbed biexponential 
sound-speed profile. 



B. Eigenfunction statistics 

Wave chaos is also refiected in eigenfunction statistics 
of the FREO. Each eigenfunction can be expressed as a 
superposition of normal modes: 



^m\^) / ^ ^mny^n \'^ } i 



(33) 



where Cmn is the m-th component of n-th eigenvector of 
the matrix G, </>„(z) is the n-th normal mode. 

There are many methods for identification of "chaotic" 
eigenfunctions. In the present work we use only one 
of them. Ray chaos can lead to intense energy ex- 
change between normal modes 0, [ij, [2J|, therefore, a 
"chaotic" eigenfunction is a compound of many modes. 
The stronger chaos, the larger number of contributing 
modes. Hence, we can characterize "chaoticity" by esti- 
mating number of principal components \4^ in the ex- 
pansion ([55]) . Number of principal components is calcu- 
lated as 



iy(n) 




(34) 



Number of principal components v is equal to 1 in an 
unperturbed waveguide and grows as scattering intesifies, 
tending asymptotically to M. 



V. ONE-STEP POINCARE MAP 

In the geometrical acoustics approximation, dynamics 
of a wavepacket is represented as motion of some ray 
bundle. The bundle remains compact in course of propa- 
gation if ray dynamics is regular. Under conditions of ray 
chaos the bundle rapidly diverges. Initial conditions for 
rays contributing to the bundle can be found by means 



of the Husimi function 



W^hb, z) = 
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(35) 



projecting an initial wavepacket $o onto phase space of 
ray equations ([7]). In the hmit fco — > oo, the transfor- 
mation (|35p turns eigenfunctions $„ of the FREO into 
some phase space sets which are invariant under the shift 
r = — >■ r = T. The procedure proposed in [22] allows 
one to find out these sets in a randomly-inhomogeneous 
waveguide. This is the one-step Poincare map (or the 
specific Poincare map). The ray analogue of the map 
looks as follows: 

K+i =p(r = T|pi,Zi), Zi+i ^q{r ^T\pi,Zi), (36) 

where p{r = T\pi,Zi) and z{r = T\pi,Zi) are the solutions 
of ray equations (O with initial conditions p{r = 0) = Pi, 
z{r = 0) = Zi. Values of p and 2, calculated at the i-th 
step of mapping, become the initial conditions for the 
[i + l)-th step. This procedure is equivalent to the usual 
Poincare map J^] for a range-periodic waveguide with 
the ray Hamiltonian 



H 



\^Uiz) 



V{z, 



(37) 



Here V{z, r) is periodic in r function 

V{z,r' + nT) = V{z,r'), < / ^ r, (38) 

where n is an integer. As it follows from (|38p . V{z^ r) 
is a sequence of identical pieces of V(z, r), each of 
them has the length r. Thus we replace the original 
randomly-perturbed Hamiltonian system by an equiv- 
alent periodically-perturbed one. This replacement is 
valid as long as we restrict ourselves by considering 
dynamics within the range interval [0 : r]. One-step 
Poincare map can be considered as the classical coun- 
terpart of the FREO. 

Owing to analogy with the usual Poincare map, the 
main property of the one-step Poincare map can be for- 
mulated as follows: each point of a continuous closed ray 
trajectory of the map fi36]) corresponds to a starting point 
of the solution of ^ which remains stable by Lyapunov 
till the range r — t. The inverse statement is not, in gen- 
eral, true. Hence, the one-step Poincare map provides a 
sufficient but not necessary criterion of stability. 

Map dMl) was studied in d [H IH^. Here we shall 
give its brief description. It is reasonable to make canon- 
ical transformation of ray variables from momentum- 
depth (p— z) to the action-angle (/,^). This procedure 
provides more appropriate representation of ray equa- 
tions. The angle variable •§ canonically conjugated to 
the action ([T0| is given by 



"^^dJ I P^""- 



(39) 



The transformed ray Hamiltonian is written as 

fl^Ho{I) + V{I,^,r). (40) 

Ray equations in terms of the new variables: 



dl _ dH _ dV 

dd dH ,^^ dV 
_ = _=c.(/) + — , 



(41) 



where w = 2tt / D is spatial frequency of a ray trajectory 
in a waveguide. Perturbation V{1, 1}) can be expanded 
into a double Fourier series 



k,k' = l 



(42) 



where il = ^-kJt. Inserting (|42|) into (|4T|) . we obtain 
AT ■ °° 



k,k' = l 



dr 2 ^^^ 01 

k.k' = l 



(43) 



dr 
If the condition 



k'D{I) = fcr 



(44) 



is fulfilled, there occurs resonance. The pair of integers 
k' and k determines multiplicity of resonance k : k' . Res- 
onances occur at certain values of the action, which cor- 
respond to the so-called resonant torus. Ray dynamics 
in a small vicinity of a resonant torus with ray action /q 
can be described using the so-called resonance approxi- 
mation [18l|, when one leaves only resonant terms in the 
r.h.s. of (|43|) . It should be mentioned that one and the 
same resonant torus corresponds to an infinite number 
of resonances with multiplicities (jfc) : (jfc'), where j is 
an integer. However, resonance Fourier amplitudes Vk,k' 
rapidly decrease with increasing k and fc', therefore, only 
few low-order resonances influence significantly ray dy- 
namics. Consequently, we can take into account only 
some finite number of dominant resonances. For further 
simplification, we make the following procedures: 

1. As y is a smooth function of z in the underwater 
sound channel considered, the derivative dV /dl is 
small compared with uj, and the sum in the second 
equation of ((43|) can be dropped out. 



2. Near resonance, spatial frequency lo can be ex- 
panded as a; = rj + WjA/, where uj'j — dw/dl. 



Then introducing new variables 

AI = I~Io, ip = ki} - k'Qr, 



(45) 



8 



and expressing Vk^k' as |Vfc^fc'| e^p{iCk,k'), we can rewrite 
(PI as 



d(A/) 



dH 



= ^lk\Vik,ik'\sm{lil; + (ik,ik') = — y-, 

1=1 ^ (46) 

where L is the number of dominant resonances, and 

H = w^^^ + Yl ^k\^ik,ik- 1 cos(ZV + Ofc.ife')- (47) 
1=1 

If L = 1, (j47|) turns into the universal Hamiltonian of 
nonhnear resonance |3J, Il8l | , and a phase space portrait of 
Eqs. (PI contains the domain of finite motion enclosed 
by the separatrix and corresponding to trapping into res- 
onance. Then maximal value A/ on the separatrix is 
determined as 



A/„ 



= 4, 



>k\Vk,. 



iOj 



(48) 



The terms with I > 1 may deform the pendulum-like 
phase space portrait and, moreover, result in the pres- 
ence of additional separatrices inside the domain of finite 
motion. The latter phenomenon can occur when the per- 
turbation oscillates with depth |22l |. 

Transition to global chaos in the one-step Poincare 
map happens when neighbouring dominant resonances 
overlap. The criterion of overlapping is the well-known 
Chirikov criterion 



A/„ 



51 



^ 1. 



(49) 



Here 51 is the distance between neighbouring dominant 
resonances in the action space. Its variability with r for 
T > D is described by equation 



51 



2-K 



(50) 



that is, increasing of r enhances resonance overlapping. 
Differences in phase space patterns corresponding to dif- 
ferent realizations of the perturbation are associated with 
phase and amplitude fluctuations of Fourier-amplitudes 
Vk^k' ■ However, the contribution of these fluctuations is 
limited, therefore, the ratio of the phase space volumes 
corresponding to regular and chaotic motion is mainly 
controlled by r and weakly varies from one realization 
to another (see below for an illustration). This property 
allows one to consider the one-step Poincare map as an 
useful tool for studying randomly-driven dynamical sys- 
tems of various physical origins. 

Resonance approximation presented above fails under 
violation of the nondegeneracy condition uj\ ^ 0. In this 
case one should use more sophisticated approaches, one 



of them is presented in [46| . We shall address this issue 
in the end of this Section. 

Fig.|4]illustrates phase space portraits calculated using 
the one-step Poincare map with three different realiza- 
tions of the sound-speed perturbation. The jump of the 
derivative dU/dz aX z — zq leads to fast growth of nu- 
merical error in ray calculations. Therefore, we replaced 
the expression ([S| for U{z) by the smoothed function 



Uiz) = Ui{z) + 



1 



1 + tanh ■ 



za 



A 



U2{z), (51) 



where A = 1 ni. Each of the phase space portraits repre- 
sents a mixed phase space structure consisted of regular 
and chaotic domains. Phase portraits with the same r 
mainly differ only in angular locations of regular islands, 
whereas their overall structure is very similar. The main 
regular domain is placed near the point z = zq, p = 
and corresponds to flat rays intersecting the horizontal 
plane with the smallest angles. This circumstance de- 
serves especial attention because stability of flat rays is 
not typical for sound propagation in the deep ocean. Nu- 
merous experiments on long-range sound propagation in 
the North-Eastern Pacific Ocean (see, for instance, [43- 
|5l|) indicate on strong irre gula rity of fiat near-axial rays, 
associated with ray chaos [52, [Sj] . The "deterministic" 
mechanism of near-axial chaos is ray scattering on ver- 
tical resonances caused by small-scale de pth oscillations 
of the sound-speed perturbation [Tj, [Tj, ^^ . These os- 
cillations are contributed from high-number modes of an 
internal- wave field. In the Sea of Japan, the effect of the 
high-number modes is weak, therefore, the sound-speed 
perturbation can be fairly described by equation (|17p . 
where depth dependence is given by a smooth function 
Yi{z). Weakness of high- number internal- wave modes is 
a peculiar feature of the Sea of Japan. It is caused by the 
specific form of the buoyancy frequency profile |28| . In 
particular, the waveguide for internal waves, determined 
by the buoyancy frequency profile, is too narrow. Con- 
sequently, it doesn't efficiently focus high-number modes 
of low-frequency internal waves. As low-frequency in- 
ternal waves give the dominant contribution into a total 
internal- wave field, depth oscillations of the perturbation 
are suppressed. 

Resonance overlapping is enhancing as r grows, and 
stable islands eventually submerge into the chaotic sea. 
However, there is a small region of stability that survives 
for distances of hundreds kilometers, transforming into a 
chain of islands around the region near z = zq, p — 0. 
This chain corresponds to the smooth minimum of the 
function D{po) depicted in Fig. [51 In the absence of inho- 
mogeneity, this minimum gives rise to the so-called shear- 
less torus in phase space, producing a weakly-divergent 
beam. It is recognized that shear less tori can possess ex- 
traordinary persistence to chaos |46l Isil - ISTt ■ therefore, a 
weakly-divergent beam can survive in the presence of a 
sound-speed perturbation. Hence, formation of a weakly- 
divergent beam can be considered as a possible mecha- 
nism responsible for unusual stability of near-axial rays. 
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observed in experiments [271, |2S 



VI. SPECTRAL STATISTICS OF FREO FOR 

THE UNDERWATER SOUND CHANNEL IN THE 

SEA OF JAPAN 



A. General remarks 



This section is devoted to numerical modeling of the 
FREO in the Sea of Jaoan. Each realization of the FREO 
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FIG. 5. Fraction of the phase space volume corresponding 
to regular motion vs distance r for various frequencies. 



was represented as a matrix in the basis of normal modes 
being solutions of the Stumi-Liouville problem (I12p with 
the boundary conditions (j22p . Only purely- water modes 
which propagate without reaching the bottom were taken 
into account. They were selected using the criterion 
Em < U{z = h), where £",„ is the m-th eigenvalue of the 
Sturm-Liouville problem (J12p . This criterion ensues from 
the WKB approximation for normal modes [9|. Number 
of trapped modes M depends on sound frequency. It is 
equal to 72 for / = 100 Hz, 179 for / = 250 Hz, 259 for 
/ = 360 Hz, and 361 for / — 500 Hz. Statistical ensem- 
bles of the FREOs, corresponding to the frequencies of 
250, 360 and 500 Hz, were calculated with 100 realiza- 
tions of the perturbation. The ensemble corresponding 
to 100 Hz was calculated with 500 realizations. For each 
realization, we constructed a family of the FREOs G{t), 
where t = 5, 10, 15, ... , 350 km. 
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FIG. 6. The Brody parameter /3 vs distance t. 



B. Eigenvalue statistics 

We calculated the ensemble-averaged level spacing dis- 
tribution p{s,t) using the formula (1321) and fitted it, for 
each value of t, by means of the Berry-Robnik (PUI) and 
Brody ([3T|) distributions. Thus, we obtained dependen- 
cies of the regular phase space volume v^- and Brody pa- 
rameter /3 on distance r. As is shown in Fig.[5l Ur rapidly 
decreases in the first 100-150 km. Then, it becomes al- 
most constant. This may indicate on the influence of the 
long-living stable islands in the vicinity of the weakly- 
divergent beam. Notably, the curves corresponding to 
250, 360 and 500 Hz are very close to each other, whereas 
the curve corresponding to 100 Hz lies above them and 
undergoes strong fluctuations which persist even with in- 
creasing number of realizations. It should be mentioned 
that the above estimates of Vr have limited accuracy be- 
cause the assumptions underlying the formula pop are 
satisfled only approximately. Therefore, the results ob- 
tained using the Berry-Robnik distribution are rather 
qualitative than quantitative, especially for low sound 
frequencies. 
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FIG. 7. Relano parameter a vs distance r. 

Approximation of level spacing statistics by means of 
the Brody distribution leads to qualitatively similar re- 
sults. They are presented in Fig. [6l For frequencies 250, 
360 and 500 Hz, the Brody parameter /? gradually grows 
from to 1, reflecting transformation from Poissonian 
to Wigner-like statistics. In the case of / = 100 Hz the 
growth is remarkably slower and affected by strong fluc- 
tuations. The fluctuations of «,■ and f3 may be induced by 
regular-to-chaotic tunneling and dynamical localization, 
the phenomena whose influence on spectral statistics is 
still not well understood. 

In addition, we used the method of spectral analysis, 
developed by A. Relano with coworkers in [53]. In this 
method, one firstly constructs a series 



Sn = X1(*^~ ^*))' 



(52) 
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where n = 1,2, . . . ,N—1, N is the total number of eigen- 
values. Then, making discrete Fourier transform 



2TTikn 



^"'= = -^E^"^-Pv A. 



one finds the power spectrum 

S{k) = 141'. 



(53) 



(54) 



Generally, ensemble-averaged spectrum obeys a power 
low 



{S{k)) ^ fc-". 



(55) 



Relano with coworkers found that regular dynamics cor- 
responds to a = 2, and global chaos results in a = 1. 
In the mixed regime, a takes on an intermediate value 
between 1 and 2 [53]. Fig. [7] demonstrates that a de- 
creases with increasing r for the frequencies of 360 and 
500 Hz, indicating gradual transition to chaos. How- 
ever, a varies relatively slowly and remains near the mid- 
dle value 1.5 for all distances considered, despite of the 
marked changes in the classical phase space portrait (see 
Fig. S]). Analogous dependencies for the frequencies of 
100 and 250 Hz exhibit strong fluctuations and, there- 
fore, are not presented in the figure. This implies that 
the method developed in |58l . |59[ provides good agree- 
ment only for relatively short wavelengths. 



C. Eigenfunction analysis 

Analysis of eigenfunctions possesses some advantages 
as compared with analysis of eigenvalues. The main ad- 
vantage is the possibility to associate each eigenfunction 
with some set of normal modes and, whereby, associate it 
with a certain geometry of propagation. We can properly 
classify eigenfunctions, taking into account the interplay 
with normal modes. Such classification can be used for 
finding wavepacket configurations whose dynamics is ex- 
pected to be less or more regular. Proper classification 
can be obtained using the parameter /x [ll| . It is defined 
as 



El 

m—l 



(56) 



In an unperturbed waveguide, only one normal mode con- 
tributes to each eigenfunction, and ^ coincides with the 
number of this mode. Taking into account the quantiza- 
tion rule (jlip . we obtain the formula 



fco 2ko 



(57) 



that gives mean action corresponding to an eigenfunc- 
tion. According to (1571) . the parameter /z determines 
phase space location of the eigenfunction and can serve 



as its identificator. Figures [8] and |9] illustrate eigenfunc- 
tion distributions in the ^-v plane, where v is number of 
principal components (p4| , for the frequencies of 500 and 
100 Hz, respectively. 

Let us firstly consider Fig. [5] corresponding to the fre- 
quency of 500 Hz. Mode coupling is relatively weak for 
small values of t, therefore, the distribution is mainly 
concentrated near v = 1. We want to focus attention on 
the dense spots elongated in the i^-direction. They cor- 
respond to eigenfunctions reflecting mode-medium res- 
onance [9|] being the wave counterpart of ray-medium 
resonance (|33]). Indeed, an isolated ray-medium reso- 
nance produces oscillations of ray action inside the in- 
terval Iq — A/max ^ I ^ !() + A/max, where /o is a res- 
onance action, A/,„ax is resonance width in the action 
space. According to the principle of ray-mode duality, 
these oscillations of action correspond to coherent tran- 
sitions between the normal modes whose numbers satisfy 
the inequality 



mo — Am ^ 771 ^ 77lo + A777, 



(58) 



where rriQ = fco/o -I- 1/2, A777 = fcoA/max- Resonance- 
induced modal transtions give rise to eigenfunctions with 
/i ~ mo and v varying from 1 to A777. As long as res- 
onance values of the action, being determined by r, are 
the same for all realizations of the perturbation, these 
eigenfunctions form vertically-elongated concentrations 
of points in the ii-v plane. Location of mode-medium 
resonances along the /Lt-axis can be found using the for- 
mula 



k'D{I= (/)) =kT, 



(59) 



where (/) is linked to /x by (|57p. Concentrations induced 
by mode-medium resonance disappear with increasing t 
due to overlapping of mode-medium resonances and de- 
localization [601 • Delocalization leads to abrupt growth 
of number of principal components. It eventually sub- 
jects all eigenfunctions in the interval between fi ~ 100 
and /i ~ 300, resulting in the "boomerang" pattern in 
the /i-7^ plane, as illustrated in Fig. |5l^d). Left and right 
ends of the "boomerang" are formed by weakly scattered 
eigenfunctions. The left end corresponds to the almost 
horizontal near-axial propagation, that is, its regularity 
can be associated with long-living stable islands in the 
one-step Poincare map. 

Eigenfunction distribution in the /i-i/ plane for the fre- 
quency of 100 Hz possesses a more complicated struc- 
ture. It is exceptionally regular for r = 10 km and 
T = 35 km, as shown in Figs, ^a) and Hj^b). Mode- 
medium resonances reveal themselves as "stalagmites". 
Each stalagmite is drawn by a family of distinct weakly 
biased lines. Contours of the most pronounced stalagmite 
for T = 35 km are somewhat disordered and smeared. 
Some traces of stalagmite-like patterns survive even for 
distances of hundreds kilometers, despite of global over- 
lapping of ray-medium resonances. Persistence of stalag- 
mites for large r indicates on the presence of eigenstates 
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localized near periodic orbits of the one-step Poincare 
map. As long as these eigenstates correspond to non- 
speading wavepackets, we can associate such localization 
with suppresion of ray chaos and recovery of regular re- 
fraction in the vicinities of the periodic orbits. An anal- 
ogous phenomenon had been reported in jl2,] for range- 
periodic waveguides. 

Besides of stalagmites, Figs. [9ja) and Mjo) illustrate 
the pattens in the form of "bridges". For instance, a 
pronounced "bridge" in the left part of Fig. |9l[a) connects 
the points fi — 5, i' — 1 and fi — 30, v ^ \. The 
eigenfunctions producing the "bridges" are consisted of 
normal modes satisfying the condition 



ko{E,n — En) 



2nl 



m > n. 



(60) 



The aforementioned "bridge" in the Fig. HJa) satisfies 
(|60l) with m = 30, n = 5, and I = 9. Condition (|60l) 
is equivalent to quantum resonance between two energy 
levels. In the ray limit 



ko{E„ 



E„ 



dE 
— > [ra — n)—— 



27r(rn — n) 
D ' 



(61) 



and condition (|60p reduces to (j44|) . Each realization 
of the FREO can yield eigenfunctions manifesting res- 
onance (IBH]) . If resonance (IBH]) corresponding to some 



numbers {I, m, n) is localized, i. e. the modes m and n 
are not affected by other resonances, then the respective 
eigenfunction is a superposition of normal modes m and 



$rcs(z) 



+ C„ 



1. 



(62) 



As this takes place, the ratio of amplitudes |cm|/|c„| is 
determined by the phase of the resonance harmonics of 
the perturbation. For the perturbation (J17p . amplitude 
of the ^-th resonance harmonics reads 



Bi 



1 

2^ 



2ir/r 



bi(r) exp 



2TTlr , , 
-i dr. 



(63) 



Phase of the resonance harmonics is a random quantity 
with uniform probability density in the range [0 : 27r]. 
Each value of the phase uniquely determines the values 
of /I and ly via the formulae 



fi = |c„ipm+ |c„pn. 



^=(|c„/ + |c„|4) 1 



(64) 



(65) 



It turnes out that quantities ^ and v are correlated for 
each eigenfunction corresponding to localized resonance 
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ferences in sound scattering for different frequencies. To 
estimate these differences quantitatively, we need a suit- 
able parameter characterizing scattering strength. For 
example, chaos-induced phase space delocalization can 
be measured by the ensemble-averaged number of princi- 
pal components, divided by the total number of trapped 
modes. As it follows from Fig. [TUl the rate of delocal- 
ization increases with increasing sound frequency. This 
implies that chaotic diffusion associated with ray chaos 
degrades with increasing sound wavelength. It can be 
thought of as a manifestation of dynamical localization 
[25| , an analogue of Anderson localization, when destruc- 
tive interference supresses wavepacket spreading. 

In practice, it is useful to know fraction of the eigen- 
function ensemble corresponding to regular propagation. 
This quantity can be regarded as an analogue of the pa- 
rameter Wr in the Berry-Robnik distribution. It can be 
estimated by means of the cumulative distribution func- 
tion 



(pO)l . This results in formation of "bridges" in the ^- 
V plane. Localization becomes violated with increasing 
T, therefore, correlation between v and ^ ceases, and 
ordered "bridges" transform into disordered clouds of 
points. 

The above analysis shows that there are qualitative dif- 



i^M= p{v')dv', 



(66) 



where p{v') is the probability density function of u. We 
can conditionally distinguish two regimes of localization: 
strong localization and moderate localization. Strong lo- 
calization implies that eigenfunction of FREO is close to 
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one of normal modes of the unperturbed waveguide. To 
select strongly-localized eigenfunctions, we can use the 
inequality v ^ 2. Dependence of F{2) on the horizon- 
tal distance r is depicted in Fig. [TT] Evidently, frac- 
tion of strongly-localized eigenfunctions is much larger 
for / = 100 Hz than for higher frequencies. The curves 
corresponding to the higher frequencies are close to each 
other. They drop down in the range of small r, then 
decreasing of F{2) becomes very slow. Slow decreasing 
can be linked to the presence of the long-living islands 
of stability in the neibourhood of the weakly-divergent 
beam. 

In the regime of moderate localization, mode coupling 
can be sufficiently strong, but an eigenfunction occupies 
relatively small phase space volume, that is, number of 



principal components is limited. We use the inequality 
V ^ O.IM as the criterion of moderate localization. As 
is demonstrated in Fig. [T^ almost all eigenfunctions cor- 
responding to the frequency of 100 Hz are moderately- 
localized. This is not the case of higher frequencies, 
when fraction of moderately-localized eigenfunctions sig- 
nificantly decreases with r. We can conclude that de- 
creasing of sound frequency can result in remarkable sup- 
pression of wavepacket diffusion in phase space. 



VII. CONCLUSION 

Spectral analysis of the finite-range evolution operator 
prompts a way to explore wave dynamics in a randomly- 
inhomogeneous waveguide by means of the quasideter- 
ministic approach, involving resonances, periodic orbits, 
phase space portraits, e.t.c. This approach was originally 
proposed in [ij, [2j]. The present work is devoted to 
its further development. We demonstrate various meth- 
ods of spectral analysis. For instance, fitting of level 
spacing statistics by means of the Berry-Robnik distri- 
bution yields approximate estimate for fraction of regu- 
larly propagating normal modes of a waveguide. Also we 
use the method developed by A. Relano with coworkers, 
and consider distribution of eigenfunctions in the ^,-v 
space. In our opinion, the latter approach is the most 
robust, albeit its analytical description is lacking. The 
most important advantage of the eigenfunction analysis 
is the possibility to study separately scattering of differ- 
ent modes of a waveguide and find out the modes cor- 
responding to regular propagation. A detailed view of 
the eigenfunction distribution in the fi-v space suggests 
that the mechanism of the chaos onset with increasing 
distance can be associated with overlapping and delocal- 
ization of mode-medium resonances. The approach based 
on the statistical analysis of level spacings by means of 
the Berry-Robnik formula gives qualitative description 
of the transition but doesn't ensure quantitative agree- 
ment. The Relano method and the method based on the 
Brody distribution also give the qualitative description 
but cannot make any quantitative estimates due to their 
semiempirical nature. 

We consider the underwater sound channel in the Sea 
of Japan as an example. Our analysis shows that al- 
most horizontal near-axial sound propagation preserves 
regularity over distances of hundreds kilometers. There 
are two factors responsible for near-axial stability. The 
first factor is the peculiar hydrological structure in the 
region considered, resulting in the absence of vertical 
oscillations of the sound-speed perturbation. This cir- 
cumstance leads to a qualitatively different scenario of 
ray chaos, as compared with the well-known acoustic ex- 
periments in the North-Eastern Pacific Ocean. The sec- 
ond factor is the formation of the weakly-divergent beam 
supported by long-living stable islands in classical phase 
space. Weakly-divergent beam occurs in the vicinity of 
the shearless torus. Non-dispersive motion of quantum 
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wavepackets near shearless tori was earlier observed in 
[61| . Also, we should emphasize that sound propagation 
with the frequency of 100 Hz reveals exceptionally high 
degree of stability, indicating on different physics of scat- 
tering for low sound frequencies. This phenomenon can 
be associated with strong dynamical localization. 

We suppose that applicability of the approach pre- 
sented in this paper is not limited to the problems of wave 
propagation in random media. An analogue of the FREO 
can be readily used for studying of noise-induced quan- 
tum transport and related phenomena. We suggest that 



combination of deterministic and statistical approaches 
should provide an insightful view on details of dynamics, 
especially in the presence of intermittency or synchro- 
nization. 
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